home *** CD-ROM | disk | FTP | other *** search
/ Space & Astronomy / Space and Astronomy (October 1993).iso / mac / programs / space / DE118I.ZIP / ASINL.C < prev    next >
C/C++ Source or Header  |  1993-01-30  |  4KB  |  229 lines
  1. /*                            asinl.c
  2.  *
  3.  *    Inverse circular sine, long double precision
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * double x, y, asinl();
  10.  *
  11.  * y = asinl( x );
  12.  *
  13.  *
  14.  *
  15.  * DESCRIPTION:
  16.  *
  17.  * Returns radian angle between -pi/2 and +pi/2 whose sine is x.
  18.  *
  19.  * A rational function of the form x + x**3 P(x**2)/Q(x**2)
  20.  * is used for |x| in the interval [0, 0.5].  If |x| > 0.5 it is
  21.  * transformed by the identity
  22.  *
  23.  *    asin(x) = pi/2 - 2 asin( sqrt( (1-x)/2 ) ).
  24.  *
  25.  *
  26.  * ACCURACY:
  27.  *
  28.  *                      Relative error:
  29.  * arithmetic   domain     # trials      peak         rms
  30.  *    IEEE     -1, 1        30000       2.7e-19     4.8e-20
  31.  *
  32.  *
  33.  * ERROR MESSAGES:
  34.  *
  35.  *   message         condition      value returned
  36.  * asin domain        |x| > 1           0.0
  37.  *
  38.  */
  39. /*                            acosl()
  40.  *
  41.  *    Inverse circular cosine, long double precision
  42.  *
  43.  *
  44.  *
  45.  * SYNOPSIS:
  46.  *
  47.  * double x, y, acosl();
  48.  *
  49.  * y = acosl( x );
  50.  *
  51.  *
  52.  *
  53.  * DESCRIPTION:
  54.  *
  55.  * Returns radian angle between -pi/2 and +pi/2 whose cosine
  56.  * is x.
  57.  *
  58.  * Analytically, acos(x) = pi/2 - asin(x).  However if |x| is
  59.  * near 1, there is cancellation error in subtracting asin(x)
  60.  * from pi/2.  Hence if x < -0.5,
  61.  *
  62.  *    acos(x) =     pi - 2.0 * asin( sqrt((1+x)/2) );
  63.  *
  64.  * or if x > +0.5,
  65.  *
  66.  *    acos(x) =     2.0 * asin(  sqrt((1-x)/2) ).
  67.  *
  68.  *
  69.  * ACCURACY:
  70.  *
  71.  *                      Relative error:
  72.  * arithmetic   domain     # trials      peak         rms
  73.  *    IEEE      -1, 1       30000       1.4e-19     3.5e-20
  74.  *
  75.  *
  76.  * ERROR MESSAGES:
  77.  *
  78.  *   message         condition      value returned
  79.  * asin domain        |x| > 1           0.0
  80.  */
  81.  
  82. /*                            asin.c    */
  83.  
  84. /*
  85. Cephes Math Library Release 2.2:  December, 1990
  86. Copyright 1984, 1990 by Stephen L. Moshier
  87. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  88. */
  89.  
  90. #ifndef NOINCS
  91. #include "mconf.h"
  92. #endif
  93.  
  94. #ifdef UNK
  95. long double asinP[] = {
  96.  3.7769340062433674871612E-3L,
  97. -6.1212919176969202969441E-1L,
  98.  5.9303993515791417710775E0L,
  99. -1.8631697621590161441592E1L,
  100.  2.3314603132141795720634E1L,
  101. -1.0087146579384916260197E1L,
  102. };
  103. long double asinQ[] = {
  104. /* 1.0000000000000000000000E0L,*/
  105. -1.5684335624873146511217E1L,
  106.  7.8702951549021104258866E1L,
  107. -1.7078401170625864261444E2L,
  108.  1.6712291455718995937376E2L,
  109. -6.0522879476309497128868E1L,
  110. };
  111. #endif
  112.  
  113. #ifdef IBMPC
  114. short asinP[] = {
  115. 0x59d1,0x3509,0x7009,0xf786,0x3ff6, XPD
  116. 0xbe97,0x93e6,0x7fab,0x9cb4,0xbffe, XPD
  117. 0x8bf5,0x6810,0xd4dc,0xbdc5,0x4001, XPD
  118. 0x9bd4,0x8d86,0xb77b,0x950d,0xc003, XPD
  119. 0x3b0f,0x9e25,0x4ea5,0xba84,0x4003, XPD
  120. 0xea38,0xc6a9,0xf3cf,0xa164,0xc002, XPD
  121. };
  122. short asinQ[] = {
  123. /*0x0000,0x0000,0x0000,0x8000,0x3fff,*/
  124. 0x1229,0x8516,0x09e9,0xfaf3,0xc002, XPD
  125. 0xb5c3,0xf36f,0xe943,0x9d67,0x4005, XPD
  126. 0xe11a,0xbe0f,0xb4fd,0xaac8,0xc006, XPD
  127. 0x4c69,0x1355,0x7754,0xa71f,0x4006, XPD
  128. 0xded7,0xa9fe,0x6db7,0xf217,0xc004, XPD
  129. };
  130. #endif
  131.  
  132. #ifdef MIEEE
  133. long asinP[] = {
  134. 0x3ff60000,0xf7867009,0x350959d1,
  135. 0xbffe0000,0x9cb47fab,0x93e6be97,
  136. 0x40010000,0xbdc5d4dc,0x68108bf5,
  137. 0xc0030000,0x950db77b,0x8d869bd4,
  138. 0x40030000,0xba844ea5,0x9e253b0f,
  139. 0xc0020000,0xa164f3cf,0xc6a9ea38,
  140. };
  141. long asinQ[] = {
  142. /*0x3fff0000,0x80000000,0x00000000,*/
  143. 0xc0020000,0xfaf309e9,0x85161229,
  144. 0x40050000,0x9d67e943,0xf36fb5c3,
  145. 0xc0060000,0xaac8b4fd,0xbe0fe11a,
  146. 0x40060000,0xa71f7754,0x13554c69,
  147. 0xc0040000,0xf2176db7,0xa9feded7,
  148. };
  149. #endif
  150.  
  151. extern long double Zero, Half, One, Two;
  152.  
  153. long double asinl(x)
  154. long double x;
  155. {
  156. long double a, p, z, zz;
  157. long double ldexpl(), sqrtl(), polevll(), p1evll();
  158. short sign, flag;
  159. extern long double PIO2L;
  160.  
  161. if( x > 0 )
  162.     {
  163.     sign = 1;
  164.     a = x;
  165.     }
  166. else
  167.     {
  168.     sign = -1;
  169.     a = -x;
  170.     }
  171.  
  172. if( a > One )
  173.     {
  174.     mtherr( "asinl", DOMAIN );
  175.     return( Zero );
  176.     }
  177.  
  178. if( a > Half )
  179.     {
  180.     zz = Half -a;
  181.     zz = ldexpl( zz + Half, -1 );
  182.     z = sqrtl( zz );
  183.     flag = 1;
  184.     }
  185. else
  186.     {
  187.     z = a;
  188.     zz = z * z;
  189.     flag = 0;
  190.     }
  191.  
  192. p = zz * polevll( zz, asinP, 5)/p1evll( zz, asinQ, 5);
  193. z = z * p + z;
  194. if( flag != 0 )
  195.     {
  196.     z = z + z;
  197.     z = PIO2L - z;
  198.     }
  199. if( sign < 0 )
  200.     z = -z;
  201. return(z);
  202. }
  203.  
  204.  
  205. extern long double PIO2L, PIL;
  206.  
  207. long double acosl(x)
  208. long double x;
  209. {
  210. long double asinl(), sqrtl();
  211.  
  212. if( x < -One )
  213.     goto domerr;
  214.  
  215. if( x < -Half) 
  216.     return( PIL - Two * asinl( sqrtl(Half*(One+x)) ) );
  217.  
  218. if( x > One )
  219.     {
  220. domerr:    mtherr( "acosl", DOMAIN );
  221.     return( Zero );
  222.     }
  223.  
  224. if( x > Half )
  225.     return( Two * asinl(  sqrtl(Half*(One-x) ) ) );
  226.  
  227. return( PIO2L - asinl(x) );
  228. }
  229.